Práctica de Aprendizaje Estadístico - Parte 1

options(warn = -1)
suppressPackageStartupMessages({
  library(plotly)
  library(dplyr)
  library(glmnet)
  library(corrplot)
  library(caret)
})

0. Carga los datos y elimina la variable TRAIN.

0.1 Configuración inicial. Carga de librerías.

Nota: Nos aseguramos de tener el paquete ElemStatLearn instalado. Si no: install.packages(lib_URL)

lib_URL = 'https://cran.r-project.org/src/contrib/Archive/ElemStatLearn/ElemStatLearn_2015.6.26.2.tar.gz'
# install.packages(lib_URL)
# install.packages("plotly")

library(plotly)
library(ElemStatLearn)
library(dplyr)

Utilizamos head() para ver las primeras filas y hacernos una idea de las variables disponibles

data("prostate", package = "ElemStatLearn")
head(prostate)

##       lcavol  lweight age      lbph svi       lcp gleason pgg45       lpsa
## 1 -0.5798185 2.769459  50 -1.386294   0 -1.386294       6     0 -0.4307829
## 2 -0.9942523 3.319626  58 -1.386294   0 -1.386294       6     0 -0.1625189
## 3 -0.5108256 2.691243  74 -1.386294   0 -1.386294       7    20 -0.1625189
## 4 -1.2039728 3.282789  58 -1.386294   0 -1.386294       6     0 -0.1625189
## 5  0.7514161 3.432373  62 -1.386294   0 -1.386294       6     0  0.3715636
## 6 -1.0498221 3.228826  50 -1.386294   0 -1.386294       6     0  0.7654678
##   train
## 1  TRUE
## 2  TRUE
## 3  TRUE
## 4  TRUE
## 5  TRUE
## 6  TRUE

attach(prostate)

0.2 Eliminación de datos TRAIN.

La columna train se utiliza en estudios previos para indicar las filas de entrenamiento.

Dado que esta práctica no la requiere, eliminamos esta columna para simplificar el conjunto de datos.

prostate = prostate %>% select(-train)
head(prostate)

##       lcavol  lweight age      lbph svi       lcp gleason pgg45       lpsa
## 1 -0.5798185 2.769459  50 -1.386294   0 -1.386294       6     0 -0.4307829
## 2 -0.9942523 3.319626  58 -1.386294   0 -1.386294       6     0 -0.1625189
## 3 -0.5108256 2.691243  74 -1.386294   0 -1.386294       7    20 -0.1625189
## 4 -1.2039728 3.282789  58 -1.386294   0 -1.386294       6     0 -0.1625189
## 5  0.7514161 3.432373  62 -1.386294   0 -1.386294       6     0  0.3715636
## 6 -1.0498221 3.228826  50 -1.386294   0 -1.386294       6     0  0.7654678

1. Exploración de los datos

1.1 ¿Cuántas variables hay?

Tras eliminar la columna train, contamos con un total de 9 variables.

1.2 ¿De qué clase son?

Utilizamos la función sapply para aplicar una función, en este caso class, a todas las columnas.

sapply(prostate, class)

##    lcavol   lweight       age      lbph       svi       lcp   gleason     pgg45 
## "numeric" "numeric" "integer" "numeric" "integer" "numeric" "integer" "integer" 
##      lpsa 
## "numeric"

Los tipos de las variables son del tipo numeric e integer, por lo que todos datos son números. Dado que son datos numéricos, podemos echar un vistazo a sus distribuciones.

par(mfrow=c(3,3))
for (column in names(prostate)) {
  hist(prostate[[column]], main=paste("Histograma de", column), xlab = column)
}

Podemos observar que los valores de la variable svi se concentran en el 0 y en el 1, indicador de que esta variable podría ser un booleano. Veamos si estos son sus únicos valores.

unique(prostate$svi)

## [1] 0 1

0 y 1 son los únicos valores presentes, por lo que confirmaríamos que se trata de una variable booleana donde 1 indicaría que el cáncer ha invadido el vesículo seminal, y 0 indicaría lo contrario.

1.3 ¿Hay una variable que correspondiente al identificador de paciente? Si es así, elimínala.

No existe ninguna variable con el identificador del paciente.

1.4 ¿Hay valores nulos en alguno de los ficheros?

Averiguamos si existe algún valor nulo en todo el conjunto de datos. En caso de que lo haya, identificamos las columnas que los contienen.

any(is.na(prostate))

## [1] FALSE

Dado que no hay ningún valor nulo en nuestro set de datos, continuamos.

1.5 ¿Están estandarizadas las variables? En este punto del análisis, ¿es necesario normalizarlas?

Tras revisar los histogramas realizados en el subapartado 2 podemos concluir que, dados los rangos de las variables y su disposición, estos no se encuentran normalizados.

1.6 ¿Por qué crees que algunas variables están es escala logarítmica?

El objetivo de poner algunas de las variables en escala logarítmica es el de facilitar su posterior análisis por varios motivos:

• Reducir la asimetría en casos en los que la distribución está muy sobrecargada en uno de sus extremos donde algunos valores excepcionales distorsionarían las estadísticas de la variable, haciendo que sus distribuciones se asemejen más a las de una distribución normal (como se observa en lcavol, lweight y lpsa).
• Hacer que algunas variables se comporten de forma lineal cuando en escala natural no lo harían.

Además, cabe añadir que en los campos de la biología y la química es común usar escalas logarítmicas cuando estas hacen más fácil la interpretación de los datos.

2. Análisis de variables categóricas

Para convertir las variables svi y gleason a categóricas usamos la función factor.

prostate$svi <- factor(prostate$svi)
prostate$gleason <- factor(prostate$gleason)

A continuación, visualizamos las distribuciones usando gráficos de barras:

par(mfrow = c(1, 2))  # Ajustar disposición gráfica
barplot(table(prostate$svi), main = "Distribución de SVI", col = "lightblue")
barplot(table(prostate$gleason), main = "Distribución de Gleason", col = "lightgreen")

Se puede observar como para la mayoría de pacientes se indica que el cáncer no ha invadido la vesícula seminal, a la par que la mayoría tiene un índice de Gleason de 6 o 7. Esto puede indicar cierta relación, ya que el índice de Gleason mide la agresividad del cáncer, por lo que índices más bajos indican cánceres menos agresivos. Podríamos aventurarnos a pensar que la mayoría de cánceres agresivos no llegarían a infectar la vesícula seminal.

En el caso de la variable Edad queremos hacer una agrupación cada 5 años, por lo que usaremos la función cut para generar esas categorías:

max.age = max(prostate$age)
min.age = min(prostate$age)
# Añadimos +5 al maximo para que el rango tome tambien el ultimo grupo
prostate$age <- cut(prostate$age, breaks = seq(from=min.age, to=max.age+5, by=5), right=FALSE)
barplot(table(prostate$age), main = "Distribución de Edad", col = "salmon", width = 3)

En el caso de la edad nos encontramos una distribución que podría recordar a la normal, centrada en los 65 años, encontrándose el grueso de la distribución entre los 61 y 70 años. Cabe añadir que encontramos casos desde los 41 hasta los 79 años.

3. Análisis de frecuencias

3.1 ¿Qué porcentaje de pacientes con la puntuación de Gleason igual a 7, presenta índice igual svi igual a 0?

# Filtramos por gleason = 7
gleason.7 <- subset(prostate, gleason == 7)
gleason.7.props <- prop.table(table(gleason.7$svi))*100

Por lo tanto, un 66% de los pacientes con índice de Gleason igual a 7 presenta un svi de 0. Podemos representar esta distribución mediante un gráfico de sectores:

fig <- plot_ly(
  labels = names(gleason.7.props),
  values = as.numeric(gleason.7.props),
  type="pie",
  textinfo = "label+percent",
  hoverinfo = "none",
  insidetextorientation = "radial",
  height = 400
) %>%
  layout(
    title = "Distribución del índice de Gleason para pacientes con svi = 0",
    showlegend = TRUE
  )
fig

3.2 ¿Qué porcentaje de pacientes con índice svi igual a 0 tiene la puntuación de Gleason igual a 7?

# Filtramos por gleason = 7
svi.0 <- subset(prostate, svi == 0)
svi.0.props <- prop.table(table(svi.0$gleason))*100

Por lo tanto, un 49% de los pacientes con un svi de 0 presenta un índice de Gleason igual a 7. Podemos visualizar la distribución del índice de Gleason para pacientes con un svi\(=0\) mediante un gráfico de sectores:

fig <- plot_ly(
  labels = names(svi.0.props),
  values = as.numeric(svi.0.props),
  type="pie",
  textinfo = "label+percent",
  hoverinfo = "none",
  insidetextorientation = "radial" ,
  height = 400
) %>%
  layout(
    title = "Distribución del índice de Gleason para pacientes con svi = 0",
    showlegend = TRUE
  )
fig

3.3 Estas dos variables, ¿son independientes?

Vistos los resultados de estos dos últimos subapartados, podemos pensar que estas variables están relacionadas. Además, ya en el apartado 2 se vio como también podría existir una posible relación entre las variables y se le trató de dar una explicación lógica.

Para probar si realmente existe esta relación entre variables, exponemos los datos a un test de independencia chi cuadrado, cuya implementación en R es simple y rápida. En este test queremos probar que las variables están relacionadas, lo cual es nuestra hipótesis \(H_1\), por lo que la hipótesis nula \(H_0\) será que estas son independientes:

contingencia.svi_gleason <- table(prostate$gleason, prostate$svi)
chi_test <- chisq.test(contingencia.svi_gleason)
chi_test

## 
##  Pearson's Chi-squared test
## 
## data:  contingencia.svi_gleason
## X-squared = 15.918, df = 3, p-value = 0.001179

Al realizar el test, obtenemos un \(p\)-valor de 0.0012. Dado que los niveles de significación más comunes son \(0.05\) y \(0.01\), y nuestro \(p\)-valor es notablemente menor, podemos afirmar con seguridad que existe una relación entre estas variables.

4. Regresión lineal simple.

Análisis de la dependencia lineal de la variable lpsa con la variable lcavol.

library(dplyr)
library(ggplot2)

Ajustamos el modelo y vemos su summary():

modelo_lineal <- lm(lpsa ~ lcavol, data = prostate)

summary(modelo_lineal)

## 
## Call:
## lm(formula = lpsa ~ lcavol, data = prostate)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.67624 -0.41648  0.09859  0.50709  1.89672 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.50730    0.12194   12.36   <2e-16 ***
## lcavol       0.71932    0.06819   10.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7875 on 95 degrees of freedom
## Multiple R-squared:  0.5394, Adjusted R-squared:  0.5346 
## F-statistic: 111.3 on 1 and 95 DF,  p-value: < 2.2e-16

4.1 Interpretación del Modelo Lineal calculado

Ecuación del Modelo:

La ecuación del modelo de regresión lineal es:

\[ \text{lpsa} = 1.50730 + 0.71932 \cdot \text{lcavol} \]

Intercepto (\(\beta_0\)):

• El intercepto de 1.50730 indica el valor estimado de lpsa cuando lcavol es igual a 0.

Pendiente (\(\beta_1\)):

• El coeficiente para lcavol es 0.71932. Esto sugiere que, por cada unidad de incremento en lcavol, se espera que lpsa aumente en aproximadamente 0.71932 unidades.

4.2 Plot de los datos junto a la recta de regresión.

Crearemos una gráfico de dispersión de lcavol vs lpsa.

ggplot(prostate, 
       aes(x = lcavol, y = lpsa)) +
  
 # Puntos de los datos
geom_point(color = "blue", 
           alpha = 0.4) +

# Recta de regresión
geom_smooth(method = "lm", 
            color = "red") + 
  
labs(title = "Relación entre lcavol y lpsa",
     x = "Log de Volumen de Cáncer (lcavol)",
     y = "Log de Antígeno Prostático Específico (lpsa)") +
  
theme_minimal()

## `geom_smooth()` using formula = 'y ~ x'

4.3 Intervalos de confianza para los coeficientes del modelo con una confianza de 0.95

Los intervalos de confianza al 95% para los coeficientes del modelo son importantes para entender la precisión de las estimaciones. A continuación se presentan los intervalos calculados:

confint(modelo_lineal, level = 0.95)

##                 2.5 %    97.5 %
## (Intercept) 1.2652222 1.7493727
## lcavol      0.5839404 0.8547004

Con esto podemos estimar que con ausencia de un tumor medido (lcavol), el psa se encontrará entre valores de 1.26 y 1.74 ng/m, en su medida logarítmica.

Ademas, podemos estimar que por cada cm (10mm) de tumor aumentado, el PSA aumenta en unidades comprendidos entre 5.8 y 8.5 ng/m.

En el análisis de regresión lineal, se ha observado que tanto el coeficiente del intercepto como el coeficiente de lcavol presentan valores \(p\) muy pequeños \((< 2e-16)\). Esto indica que podemos rechazar la hipótesis nula \((H_0)\) de que los “verdaderos” coeficientes son nulos. Específicamente, esto significa que existe una relación estadísticamente significativa entre el volumen del cáncer (lcavol) y el nivel de antígeno específico de la próstata (lpsa).

4.4 RSE y la eficacia del modelo.

El Error Estándar de los Residuos (RSE) es una medida de la cantidad de variabilidad de los residuos (la diferencia entre los valores observados y los valores predichos por el modelo) en un modelo de regresión lineal. Se utiliza para evaluar la precisión de las predicciones del modelo. Un RSE más bajo indica que el modelo tiene un mejor ajuste a los datos.

Matemáticamente, el RSE se define como:

\[ RSE = \sqrt{\frac{1}{n-2} \sum (y_i - \hat{y}_i)^2} \]

donde:

• \(y_i\) son los valores observados,

• \(\hat{y}_i\) son los valores predichos por el modelo,

• \(n\) es el número de observaciones.

Si echamos un vistazo de nuevo a nuestro resultado:

summary(modelo_lineal)

## 
## Call:
## lm(formula = lpsa ~ lcavol, data = prostate)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.67624 -0.41648  0.09859  0.50709  1.89672 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.50730    0.12194   12.36   <2e-16 ***
## lcavol       0.71932    0.06819   10.55   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7875 on 95 degrees of freedom
## Multiple R-squared:  0.5394, Adjusted R-squared:  0.5346 
## F-statistic: 111.3 on 1 and 95 DF,  p-value: < 2.2e-16

RSE

Un RSE de 0.7875 indica que, en promedio, las predicciones del modelo se desvían de los valores observados de la variable de respuesta (lpsa) en aproximadamente 0.7875 unidades en la escala logarítmica.

Indicando que por cada cm(10mm) de aumento en el lcavol tenemos un desvía de 7.8 ng/m.

R-cuadrado (R²)

El coeficiente de determinación \(R^2\) es una métrica que indica la proporción de la varianza en la variable dependiente lpsa que es predecible a partir de la variable independiente lcavol. En este modelo, el \(R^2\) es:

\(R^2 = 0.5394\)

• Interpretación: Aproximadamente el 53.94% de la variabilidad en los niveles de PSA puede ser explicada por el volumen de cáncer. Esto sugiere una relación moderada entre las variables.

R-cuadrado Ajustado (Adjusted R²)

El \(R^2\) ajustado, que en este caso es:

\(R^2\) Ajustado = 0.5346

• Interpretación: El \(R^2\) ajustado penaliza el \(R^2\) por el número de predictores en el modelo, lo que lo hace más adecuado cuando se comparan modelos con diferentes números de variables. Este valor indica que el modelo tiene un ajuste similar al \(R^2\), sugiriendo que el volumen de cáncer sigue siendo un predictor significativo.

4. 6 Interpretación del Modelo Lineal (psa y cavol)

Ecuación del Modelo

El modelo de regresión lineal simple que hemos calculado es:

\[ \text{lpsa} = 1.50730 + 0.71932 \times \text{lcavol} \]

Relación entre PSA y cavol

1. Modelo Original: La ecuación del modelo de regresión lineal que obtuvimos es:

\[ \text{lpsa} = 1.50730 + 0.71932 \times \text{lcavol} \]

\[ ln(\text{psa}) = 1.50730 + 0.71932 \times ln(\text{cavol}) \]

2. Despejar PSA: Como \(\text{lpsa}\) es el logaritmo natural del PSA, puedes deshacerte del logaritmo aplicando la función exponencial:

\[ \text{PSA} = e^{ln(\text{psa})} = e^{(1.50730 + 0.71932 \times ln(\text{cavol}))} \]

\[ \text{PSA} = e^{1.50730 } + e^{(0.71932 \times ln(\text{cavol}))} \]

\[ \text{PSA} = e^{1.50730 } + (e^{ln(\text{cavol})})^{0.71932} \]

3. Transformación a cavol: Para expresar esto en términos de cavol (el volumen de cáncer en su forma original), tienes que revertir el logaritmo:

\[ \text{cavol} = e^{ln(\text{cavol})} \]

Ahora sustituimos lcavol por su expresión en términos de cavol:

\[ \text{PSA} = e^{1.50730 } +\text{cavol}^{0.71932} \]

4. Relación Final: La relación entre PSA y cavol queda como:

\[ \text{PSA} = C \cdot \text{cavol}^{b} \] Donde \(C = e^{1.50730}\) y \(b = 0.71932\). Esto significa que PSA es proporcional a cavol elevado a un exponente.

Interpretación de la Relación

• C: Este es un coeficiente constante que representa el nivel de PSA cuando cavol es igual a 1 (o la unidad de referencia).
• b (0.71932): Indica que por cada aumento en el volumen de cáncer (cavol), el PSA aumenta en un porcentaje específico que depende del valor de cavol.

5. Regresión lineal múltiple.

5.1 Observación de correlaciones entre variables.

library(corrplot)
prostate$svi = as.numeric(prostate$svi)

prostate2 = prostate

prostate2$age_num <- as.numeric(prostate2$age)
prostate2$pgg45_num <- as.numeric(prostate2$pgg45)
prostate2$gleason_num <- as.numeric(prostate2$gleason)

cor_data = prostate2 %>% select(-age, -pgg45, -gleason)
cor_matrix = cor(cor_data)

corrplot(cor_matrix, 
         diag = FALSE,
         method = "circle", 
         type = "upper", 
         tl.col = "black", 
         tl.srt = 25,
         addCoef.col = "black")

La matriz de correlación muestra las correlaciones de Pearson entre las variables del conjunto de datos, con valores que van de -1 a 1.

Análisis de la Matriz de Correlación

Se ha realizado un análisis de correlación entre las variables del conjunto de datos. A continuación se presentan los resultados y las interpretaciones:

1. (lcavol, lpsa)
   • Máxima correlación mutua: (0.734).
   • Interpretación: La fuerte correlación entre el volumen del cáncer en escala logarítmica (lcavol) y el nivel de antígeno específico de la próstata (PSA) sugiere que a medida que aumenta el volumen del cáncer, también tienden a aumentar los niveles de PSA.
2. (lbph, lweight):
   • Máxima correlación mutua: (0.442).
   • Interpretación: La correlación entre el peso en escala logarítmica (lweight) y la hipertrofia benigna de próstata (lbph) indica que hay una relación positiva entre ambos. Esto sugiere que los hombres con mayor peso pueden tener mayor riesgo de desarrollar condiciones benignas en la próstata.
3. lcp:
   • Máxima correlación: lcavol (0.675).
   • Interpretación: La correlación entre el nivel de extensión del cáncer dentro de la cápsula (lcp) y el volumen del cáncer (lcavol) tienen trivial relación.
4. svi:
   • Máxima correlación: lcp (0.68).
   • Interpretación: La correlación entre el nivel de invasión seminal (svi) y el nivel de extensión del cáncer dentro de la cápsula (lcp) es positiva y moderada (0.68). Esto indica que a medida que aumenta la invasión a las vesículas seminales, también podría aumentar la extensión del cáncer dentro de la cápsula.
5. age_num:
   • Máxima correlación: lbph (0.364).
   • Interpretación: La correlación moderada entre la edad en formato numérico (age_num) y la hipertrofia prostática benigna (lbph) sugiere que a mayor edad, la probabilidad de presentar hipertrofia benigna aumenta.
6. (pgg45_num, gleason_num):
   • Máxima correlación mutua: (0.752).
   • Interpretación: La alta correlación entre el porcentaje de patrón Gleason 4 o 5 (pgg45_num) y la puntuación total de Gleason (gleason_num) refuerza que un mayor porcentaje de patrones más agresivos se asocia con una puntuación de Gleason más alta.

5.2 Modelo de regresión lineal para predecir lpsa.

modelo_lpsa = lm(lpsa ~ lcavol + lweight + lbph + lcp + svi, 
                  data = prostate)

summary(modelo_lpsa)

## 
## Call:
## lm(formula = lpsa ~ lcavol + lweight + lbph + lcp + svi, data = prostate)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.84878 -0.38372 -0.00413  0.45189  1.55468 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.12471    0.73259  -1.535  0.12819    
## lcavol       0.54790    0.08637   6.343 8.56e-09 ***
## lweight      0.53036    0.19769   2.683  0.00867 ** 
## lbph         0.07999    0.05643   1.418  0.15971    
## lcp         -0.03638    0.08088  -0.450  0.65391    
## svi          0.75975    0.24122   3.150  0.00221 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7071 on 91 degrees of freedom
## Multiple R-squared:  0.6443, Adjusted R-squared:  0.6248 
## F-statistic: 32.97 on 5 and 91 DF,  p-value: < 2.2e-16

Siendo la siguiente su función:

\[ \text{lpsa} = -1.12471 + 0.54790 \cdot \text{lcavol} + 0.53036 \cdot \text{lweight} + 0.07999 \cdot \text{lbph} - 0.03638 \cdot \text{lcp} + 0.75975 \cdot \text{svi} + \epsilon \]

Coeficientes:

• lcavol: 0.54790 (muy significativo, p = 8.56e-09). Esto indica que, manteniendo las otras variables constantes, un aumento de una unidad en el volumen del cáncer (lcavol) se asocia con un aumento de 0.54790 unidades en los niveles de PSA (lpsa).

• lweight: 0.53036 (significativo, p = 0.00867). Similar a lcavol, un aumento de una unidad en el peso del paciente (lweight) se asocia con un aumento en los niveles de PSA, sugiriendo que el peso del paciente podría influir en el diagnóstico.

• lbph: 0.07999 (no significativo, p = 0.15971). La hipertrofia benigna de próstata (lbph) no tiene un efecto significativo en los niveles de PSA, ya que su valor de p es mayor que 0.05.

• lcp: -0.03638 (no significativo, p = 0.65391). El nivel de extensión del cáncer dentro de la cápsula (lcp) no tiene un efecto significativo sobre los niveles de PSA.

• svi: 0.75975 (muy significativo, p = 0.00221). La invasión seminal (svi) tiene una asociación significativa con los niveles de PSA, indicando que a medida que la invasión seminal aumenta, también aumentan los niveles de PSA.

Rendimiento

\(R^{2}\) ajustado: 0.6248. El modelo explica aproximadamente el 62.48% de la variación en los niveles de PSA. Esto sugiere una mejora en la capacidad predictiva con la adición de la variable svi, aunque aún hay un 37.52% de la variación que no se explica por las variables incluidas en el modelo.

RSE (Error estándar de los residuos): 0.7071. Este valor indica el tamaño promedio de los errores en las predicciones del modelo, siendo relativamente bajo, lo que sugiere que las predicciones son bastante precisas.

F-estadístico 32.97 con p-valor < 2.2e-16. Este resultado indica que el modelo es altamente significativo en su conjunto, lo que significa que al menos una de las variables predictoras tiene un impacto significativo en los niveles de PSA.

Conclusiones

El modelo de regresión muestra que lcavol, lweight, y svi son los factores más relevantes para predecir los niveles de PSA, con relaciones significativas y positivas. En cambio, lbph y lcp no tienen un impacto significativo en los niveles de PSA. El \(R^{2}\) ajustado de 0.6248 indica que el modelo explica el 62.48% de la variabilidad en PSA, y el RSE bajo sugiere predicciones precisas. El modelo es útil, pero podría mejorarse considerando otro conjunto de variables.

6. Modelo de Ridge y Lasso

Para realizar las predicciones usando los modelos de Ridge y Lasso utilizaremos la librería glmnet. Primero preparamos nuestro conjunto de datos lo dividimos entre el conjunto de training y el de testeo. Para crear nuestro training set utilizaremos el 70% de las entradas de nuestro conjunto de datos.

library(glmnet)

# Preparamos un dataset sin las variables que no usamos
datos_usados <- subset(prostate, select = -c(age, pgg45, gleason))

# Dividimos el dataset entre predictores y el valor a predecir
X <- as.matrix(datos_usados[, -which(names(datos_usados) == "lpsa")])
y <- datos_usados$lpsa

# Creamos los conjuntos de train y test, especificando la seed para poder reproducir la ejecución
# con los mismos sets
set.seed(123) 
train_idx <- sample(1:nrow(X), size = 0.7 * nrow(X))
X.train <- X[train_idx,]
y.train <- y[train_idx]
X.test <- X[-train_idx,]
y.test <- y[-train_idx]

6.1 Modelo de Ridge

Para utilizar el modelo de Ridge con glmnet debemos poner el parámetro \(\alpha = 0\). Utilizaremos un grid con 200 valores de \(\lambda\) diferentes, desde \(10^{4}\) hasta \(10^{-4}\). Estos valores se han determinado después de realizar la gráfica varias veces para asegurarnos de que las regiones que en ella se observan tienen una representación similar. A partir del modelo, podemos obtener el conjunto de coeficientes y los lambdas correspondientes a cada uno de los 200 modelos obtenidos.

lambdas = 10^seq(4,-4, length = 200)

ridge.model <- glmnet(X.train, y.train, alpha = 0, lambda = lambdas)
# Obtenemos los coeficientes y los lambdas
ridge.coefs <- as.matrix(ridge.model$beta)
ridge.lambda <- ridge.model$lambda

matplot(log(ridge.lambda), t(ridge.coefs), type = "l", lty = 1, col = rainbow(nrow(ridge.coefs)),
        xlab = "log(λ)", ylab = "Coeficientes", main = "Ridge: Coeficientes vs log(λ)")
legend("topright", legend = colnames(X), col = rainbow(nrow(ridge.coefs)), lty = 1, cex = 0.6)

Vemos claramente 3 regiones en la gráfica. Una primera región a la izquierda, hasta \(\log\lambda \approx -4\), donde los coeficientes son constantes. Una segunda región intermedia, aproximadamente entre \(-4\) y \(4\) donde estos valores varían, disminuyendo todos los coeficientes a excepción de lcp, el cual comienza con un valor negativo y crece hasta tener uno positivo. Finaliza con una tercera región, donde los coeficientes convergen a un valor cercano a 0.

Obtenemos el valor óptimo estimado de \(\lambda\) usando validación cruzada, con la función cv.glmnet con el parámetro \(\alpha=0\):

# La seleccion hecha en el CV es tambien aleatoria, asi que debebemos
# especificar la seed para asegurarnos de que podemos reprocir la prueba
set.seed(123) 
# Obtenemos el mejor valor de lambda usando cross-validation
cv.ridge.model <- cv.glmnet(X.train, y.train, alpha=0, lambda=lambdas)
lambda.min.ridge <- cv.ridge.model$lambda.min

• Valor óptimo: \(\lambda_\text{Ridge} = 0.2866\).
• Valor óptimo: \(\log\lambda_\text{Ridge} = -1.2496\).

El valor óptimo del \(\log\lambda_\text{Ridge}\) nos permite ubicar el valor del \(\lambda\) óptimo en la gráfica anterior, ubicándola en la región intermedia descrita con anterioridad.

A continuación presentamos la gráfica del \(MSE\) frente a \(\log\lambda\) correspondiente al modelo de Ridge.

plot(cv.ridge.model)

Vamos ahora a hacer la predicción de lpsa usando el conjunto de testeo.

pred.ridge <- predict(cv.ridge.model, s = lambda.min.ridge, newx = X.test)

6.2 Modelo de Lasso

Para utilizar el modelo de Lasso con glmnet debemos poner el parámetro \(\alpha = 1\). Utilizaremos el mismo grid de usado para el modelo de Ridge, pero con valores de \(\lambda\) entre \(10^1\) y \(10^{-4}\), valores que obtenemos de igual manera tras reproducir la gráfica varias veces. Después procedemos de forma similar.

lambdas = 10^seq(1,-4, length = 200)
lasso.model <- glmnet(X.train, y.train, alpha = 1, lambda = lambdas)
# Obtenemos los coeficientes y los lambdas
lasso.coefs <- as.matrix(lasso.model$beta)
lasso.lambda <- lasso.model$lambda

matplot(log(lasso.lambda), t(lasso.coefs), type = "l", lty = 1, col = rainbow(nrow(lasso.coefs)),
        xlab = "log(λ)", ylab = "Coeficientes", main = "Lasso: Coeficientes vs log(λ)")
legend("topright", legend = colnames(X), col = rainbow(nrow(lasso.coefs)), lty = 1, cex = 0.6)

Podemos observar un gráfico similar al obtenido para Ridge, con la diferencia de que Lasso sí elimina predictores, es decir hace sus coeficientes 0. Contando esta diferencia, vemos también una representación análoga de las 3 regiones que vimos en el gráfico del modelo de Ridge, con una primera región donde el coeficiente de lcp es negativo, hasta que se elimina, una región intermedia donde se regulan el resto de parámetros, y una final donde se eliminan todos los coeficientes (en Ridge también tendían a 0). Es por esto que podemos suponer que, también en este caso, el valor óptimo de \(\lambda\) se encontrará en la región intermedia.

También podemos visualizar el plot por defecto del modelo de Lasso, el cual nos muestra la evolución de los coeficientes con el L1 Norm, que representa el error que se intenta minimizar al utilizar este modelo:

\[ L1_\text{Norm} = RSS + \lambda \sum ^p_{j=1}|\beta_j| \]

plot(lasso.model, col = rainbow(nrow(lasso.coefs)))
legend("topright", legend = colnames(X), col = rainbow(nrow(lasso.coefs),start=0), lty = 1, cex = 0.6)

Obtenemos el valor óptimo estimado de \(\lambda\) usando validación cruzada, con la función cv.glmnet con el parámetro \(\alpha=1\):

# La seleccion hecha en el CV es tambien aleatoria, asi que debebemos especificar la 
# seed para asegurarnos de que podemos reprocir la prueba
set.seed(123)
# Obtenemos el mejor valor de lambda usando cross-validation
cv.lasso.model <- cv.glmnet(X.train, y.train, alpha=1, lambda=lambdas)
lambda.min.lasso <- cv.lasso.model$lambda.min
log(lambda.min.lasso)

## [1] -2.267873

• Valor óptimo: \(\lambda_\text{Ridge} = 0.1035\).
• Valor óptimo: \(\log\lambda_\text{Ridge} = -2.2679\).

A continuación presentamos la gráfica del \(MSE\) frente a \(\log\lambda\) correspondiente al modelo de Lasso.

plot(cv.lasso.model)

Vamos ahora a hacer la predicción de lpsa usando el conjunto de testeo.

pred.lasso <- predict(cv.lasso.model, s = lambda.min.lasso, newx = X.test)

6.3 Comparación con modelo lineal múltiple

1. Coeficientes: Comparamos los coeficientes obtenidos en cada modelo.

Creamos un data.frame que contenga los coeficientes para cada modelo. Esta tabla nos sirve tanto para poder comparar cuantitativamente los valores de cada coeficiente, como para preparar los datos para presentarlos en un gráfico usando ggplot.

coef.multi <- coef(modelo_lpsa)[-1]
coef.ridge <- coef(cv.ridge.model, s=lambda.min.ridge)[-1]
coef.lasso <- coef(cv.lasso.model, s=lambda.min.lasso)[-1]

coef.data <- data.frame(
  variable = names(coef.multi),
  multiple = as.numeric(coef.multi),
  ridge = as.numeric(coef.ridge),
  lasso = as.numeric(coef.lasso)
)

coef.data

##   variable    multiple      ridge     lasso
## 1   lcavol  0.54789698 0.33723614 0.3963463
## 2  lweight  0.53035882 0.55018591 0.4842218
## 3     lbph  0.07999311 0.04808135 0.0000000
## 4      lcp -0.03638017 0.57440441 0.5478204
## 5      svi  0.75974935 0.06876884 0.0000000

Cambiamos la tabla de datos a formato wide para poder mapearlo correctamente usando ggplot

coef.data.long <- reshape2::melt(coef.data, id.vars = "variable", 
                                 variable.name = "Modelo", value.name = "Coeficiente")

ggplot(coef.data.long, aes(x = variable, y = Coeficiente, fill = Modelo)) +
  geom_bar(stat = "identity", position = "dodge") + 
  labs(title = "Comparación de Coeficientes entre Modelos de Regresión",
    x = "Predictor", y = "Peso del Coeficiente") +
  theme_minimal()

• lbph: Lo primero que salta a la vista es que el modelo Lasso ha eliminado este predictor a diferencia del resto, ya que el modelo Lasso es el único que realiza selección de variables. De forma similar, el resto de modelos le otorgan un peso muy bajo.

• lcavol: Los tres modelos coinciden en otorgar un gran peso a este predictor, destacando el modelo lineal. Lasso y Ridge otorgan un valor similar entre ellos, siendo el peso en el modelo de Lasso mayor.

• lcp: Este caso es especialmente llamativo, y es que mientras los modelos Lasso y Ridge coinciden en otorgarle un peso alto (es el coeficiente de mayor valor en ambos), en el modelo de regresión múltiple este valor no es solo mucho menor, sino que además es negativo.

• lweight: En este predictor encontramos el mayor consenso entre los tres modelos, otorgándole un gran peso los tres modelos.

• svi: Este es otro caso de gran contraste entre los modelos, y es que mientras es el valor con mayor peso en el modelo múltiple, con una gran distancia al segundo, el modelo de Ridge le otorga un valor mínimo y el modelo de Lasso directamente lo elimina.

El comportamiento tan diferente de los coeficientes lcp y svi podría indicar la existencia de una correlación de importancia notable entre dichas variables.

2. Rendimiento:

Vamos a calcular la predicción del modelo lineal (el resto ya han sido calculados), y a continuación calcularemos y compararemos el MSE (error cuadrático medio) que resulta de cada modelo.

# Funciones para calcular las métricas
mse <- function(y,y2) {
  se = (y-y2)^2
  mean(se)
}

rsq <- function(y,y2) {
  rss = sum((y-y2)^2)
  tss = sum((mean(y)-y)^2)
  1 - rss / tss
}

rsq.adj <- function(y,y2,p) {
  n = length(y)
  1 - ((1-rsq(y,y2))*(n-1))/(n-p-1)
}

pred.multiple <- predict(modelo_lpsa, newdata = as.data.frame(X.test))

resultados <- data.frame(
  Metric = c("MSE", "R2", "R2adj"),
  Multiple = c(
    mse(y.test, pred.multiple), 
    rsq(y.test, pred.multiple), 
    rsq.adj(y.test, pred.multiple, 5)
  ),
  Ridge = c(
    mse(y.test, pred.ridge), 
    rsq(y.test, pred.ridge), 
    rsq.adj(y.test, pred.ridge, 5)
  ), 
  Lasso = c(
    mse(y.test, pred.lasso),
    rsq(y.test, pred.lasso),
    rsq.adj(y.test, pred.lasso, 4) # p = 4 aquí ya que, como recordamos, lasso elimina lbph
  )
)

resultados

##   Metric  Multiple     Ridge     Lasso
## 1    MSE 0.4894444 0.6175785 0.6350265
## 2     R2 0.7134772 0.6384670 0.6282529
## 3  R2adj 0.6537850 0.5631476 0.5687733

Replicamos lo hecho para los coeficientes para graficar los resultados obtenidos:

resultados.long <- reshape2::melt(resultados, id.vars = "Metric", variable.name = "Modelo", value.name = "Valor")

ggplot(resultados.long, aes(Metric, Valor, fill = Modelo)) +
  geom_bar(stat= "identity", position="dodge") +
  labs(title = "Comparación de Métricas de Rendimiento entre Modelos",
    x = "Métrica", y = "Valor de la Métrica") +
  theme_minimal()

6.4 Conclusiones

Como podemos observar, el modelo que mejor ha funcionado en nuestro caso es el de la regresión múltiple, teniendo un resultado mejor que el de los otros modelos tras medir su rendimiento usando MSE, \(R^2\) y \(R^2_\text{adj}\) en el conjunto de testeo.

Así, podemos concluir que los modelos que penalizan un número mayor de predictores, en este caso, funciona peor que el modelo de regresión múltiple que no cuenta con esa penalización. Además, tampoco representan una mejora en el rendimiento o en el espacio de almacenamiento, dado que contamos con un número muy limitado de predictores y de observaciones, por lo que no se recomendaría su uso en situaciones similares.

Para responder a si los modelos lineales podrían ser utilizados o no, revisamos los valores de las métricas obtenidas para el modelo de regresión múltiple:

• \(R^2 = 0.7134772\)
• \(R^2_\text{adj} = 0.653785\)

Como vemos, nuestro modelo funciona relativamente bien con el conjunto de testeo, por lo que podemos concluir que este modelo podría llegar a ser utilizado para realizar predicciones.

7. LDA

suppressPackageStartupMessages({
  library(MASS)
})

Ajustar el modelo LDA utilizando svi como variable de respuesta y lcavol, lcp, lpsa como variables predictoras.

prostate$svi <- as.factor(prostate$svi)
lda_model <- lda(svi ~ lcavol + lcp + lpsa, data = prostate)

7.1 Visualizar los Coeficientes del Modelo

• El modelo muestra las probabilidades de cada grupo (0 y 1). Estas indican la proporción de datos que pertenece a cada clase.

• Las medias de grupo presentan los valores medios de las variables predictoras (lcavol, lcp, lpsa) para cada clase del objetivo svi.

• Los coeficientes de los discriminantes lineales reflejan la importancia de cada variable en la construcción del discriminante lineal.

lda_model

## Call:
## lda(svi ~ lcavol + lcp + lpsa, data = prostate)
## 
## Prior probabilities of groups:
##         1         2 
## 0.7835052 0.2164948 
## 
## Group means:
##     lcavol        lcp     lpsa
## 1 1.017892 -0.6715458 2.136592
## 2 2.551959  1.6018579 3.715360
## 
## Coefficients of linear discriminants:
##                LD1
## lcavol -0.08659598
## lcp     0.76640382
## lpsa    0.53153340

7.2 Predicciones con el Modelo LDA

• class: contiene las clases predichas por el modelo.

• posterior: contiene las probabilidades posteriores de cada clase para cada observación.

• x: valores de los discriminantes lineales para cada observación.

lda_model.pred <- predict(lda_model)
names(lda_model.pred)

## [1] "class"     "posterior" "x"

7.3 Observar las Clases Predichas

Las primeras tres observaciones fueron clasificadas en la clase 0, lo que indica que el modelo ha predicho la ausencia de svi para estos casos.

lda_model.pred$class[1:3]

## [1] 1 1 1
## Levels: 1 2

7.4 Matriz de Confusión

• La matriz de confusión muestra que el modelo clasificó correctamente 69 observaciones del grupo 0 y 17 del grupo 1.

• Hay 7 casos del grupo 1 que fueron clasificados incorrectamente como 0, y 4 casos del grupo 0 clasificados incorrectamente como 1.

Esto podría estar relacionado con un desbalance en las clases, ya que el grupo 0 tiene una representación mucho mayor (aproximadamente el 78%) que el grupo 1, lo que puede hacer que el modelo esté más inclinado a predecir el grupo mayoritario.

Este desbalance puede afectar la capacidad del modelo para detectar correctamente los casos del grupo minoritario, lo que es crucial en escenarios donde los casos de la clase menos frecuente son más importantes, como en la detección de enfermedades.

library(plotly)

conf_matrix <- table(lda_model.pred$class, svi)

conf_matrix_df <- as.data.frame(as.table(conf_matrix))

heatmap <- plot_ly(
  data = conf_matrix_df,
  x = ~Var1,
  y = ~svi,
  z = ~Freq,
  type = "heatmap",
  colors = c("white", "blue"),
  showscale = FALSE
) %>%
  layout(
    title = "Predicciones de SVI",
    xaxis = list(title = "Predicción"),
    yaxis = list(title = "Clase Real")
  )

heatmap %>%
  add_annotations(
    x = conf_matrix_df$Var1,
    y = conf_matrix_df$svi,
    text = conf_matrix_df$Freq, 
    showarrow = FALSE,
    font = list(size = 12, color = "black")
  )

7.5 Probabilidades Posteriores.

• Las probabilidades posteriores indican cuán probable es que cada observación pertenezca a cada una de las clases (0 o 1).

• En este caso, para las primeras tres observaciones, la probabilidad de pertenecer a la clase 0 es extremadamente alta (cerca de 1).

lda_model.pred$posterior[1:3, ]

##           1            2
## 1 0.9998211 0.0001789316
## 2 0.9997230 0.0002769794
## 3 0.9997500 0.0002500003

7.6 Valores de los Discriminantes Lineales

• Los valores de los discriminantes lineales representan las puntuaciones de cada observación en la función discriminante.

• Estos valores son útiles para entender cómo el modelo discrimina entre los grupos basándose en las variables predictoras.

lda_model.pred$x[1:3]

## [1] -2.304200 -2.125721 -2.167584

7.7 Conclusiones

El modelo LDA ha mostrado un rendimiento desigual debido a un posible desbalanceo de clases, donde el grupo 0 es significativamente más grande que el grupo 1. Esto ha llevado a que el modelo clasifique correctamente 69 observaciones del grupo 0, pero solo 17 del grupo 1. Aunque el número de falsos positivos es relativamente bajo (4), los falsos negativos (7) indican que el modelo no está detectando adecuadamente todas las observaciones del grupo 1. Este comportamiento es común en escenarios con clases desbalanceadas.

Aunque el modelo no es completamente ineficaz, su capacidad para identificar el grupo minoritario es limitada. Se podrían explorar técnicas como el ajuste el balanceo de clases (por ejemplo, sobremuestreo o submuestreo) para mejorar el rendimiento en la clasificación del grupo 1 y obtener una predicción más equilibrada.

8. Regresión Logística

8.1 svi en función de lcavol, lcp y lpsa

Veamos primero las distribuciones del svi frente a los predictores propuestos:

boxplot1 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lcavol, 
                    type = "box", name="lcavol", width = 900) %>%
  layout(xaxis = list(title = "svi"), yaxis = list(title = "lcavol"))

boxplot2 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lcp, 
                    type = "box", name="lcp", width = 900) %>%
  layout(xaxis = list(title = "svi"), yaxis = list(title = "lcp"))

boxplot3 <- plot_ly(data = prostate, x = ~factor(svi), y = ~lpsa, 
                    type = "box", name="lpsa", width = 900) %>%
  layout(xaxis = list(title = "svi"), yaxis = list(title = "lpsa"))

fig <- subplot(boxplot1, boxplot2, boxplot3, nrows = 1, titleX = TRUE, titleY = TRUE) %>%
  layout(title="Distribución de los predictores según svi")
fig

Como podemos observar, la distribución de svi respecto al resto de variables es muy diferenciable. A simple vista podemos ver una clara relación entre valores negativos del svi con los tres predictores propuestos.

Vamos ahora a calcular el modelo. Dado que la variable svi es binaria (1 o 0), especificamos el parámetro family=binomial:

# Nos aseguramos de que la variable svi esté en formato factor, por si se ha cambiado en algun apartado anterior
prostate$svi <- as.factor(prostate$svi)

modelo_logistico <- glm(svi ~ lcavol + lcp + lpsa, data = prostate, family = binomial)
# Creamos una lista con los coeficientes preparados para el display
coefs = sapply(modelo_logistico$coefficients, function(x) round(x,2))

El modelo resulta en la siguiente función:

\[ \text{svi} = -8.58 + 0.12 \cdot \text{lcavol} + 1.31 \cdot \text{lcp}+ 2.14 \cdot \text{lpsa} \]

8.2 Conclusiones sobre el modelo

Veamos el summary del modelo:

summary(modelo_logistico)

## 
## Call:
## glm(formula = svi ~ lcavol + lcp + lpsa, family = binomial, data = prostate)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -8.5845     2.4942  -3.442 0.000578 ***
## lcavol        0.1209     0.7146   0.169 0.865645    
## lcp           1.3113     0.4240   3.093 0.001984 ** 
## lpsa          2.1431     0.8037   2.667 0.007662 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 101.353  on 96  degrees of freedom
## Residual deviance:  39.624  on 93  degrees of freedom
## AIC: 47.624
## 
## Number of Fisher Scoring iterations: 7

Lo primero que salta a la vista es que el error estimado para la variable lcavol es unas 6 veces mayor que el coeficiente en sí, lo cual nos indica que podría ser que la significancia de este predictor sea baja. Además, recordemos del apartado 5 que lcavol guardaba una alta correlación con las otras variables, lcp y lpsa. Esto se nos confirma al observar el \(p-\)valor de esta variable, de aproximadamente \(0.87\), lo cual nos indica que este predictor no es estadísticamente significativo.

Vamos a realizar una predicción de nuestro set original usando el modelo:

threshold = 0.5
log.probs = predict(modelo_logistico, type = "response")
log.preds = rep(0, length(log.probs))
log.preds[log.probs > threshold] = 1

Veamos la tabla de valores originales:

table(log.preds)

## log.preds
##  0  1 
## 79 18

Y a continuación la tabla de valores obtenidos:

table(prostate$svi)

## 
##  1  2 
## 76 21

Como se puede observar, el modelo solo falló en 3 ocasiones, dando un \(0\) cuando debería ser \(1\), lo cual indica que el modelo hizo un buen trabajo al predecir el conjunto de entrenamiento, con un 97% de precisión. Es importante hacer notar que el dataset está cargado hacia svi = 0, es decir, hay muchas más observaciones en las que svi = 1, lo cual en general no es bueno para la generación de modelos, y en particular puede ser que haya influenciado a esta tendencia al 0 que muestra nuestro modelo.

De cualquier forma, en esta ocasión no hemos dividido el dataset en training y testing, por lo que estamos examinando el rendimiento de nuestro modelo en el mismo conjunto de datos con el que lo hemos entrenado. En este caso habría que tener en cuenta dos cosas:

1. La precisión del modelo no es representativa ya que, como decimos, no estamos utilizando un conjunto de testeo.
2. Al analizar la precisión y obtener un resultado tan bueno, deberíamos tener cuidado con el overfitting.

Por estos motivos, se concluye que, a priori, no se puede decir que sea un buen modelo sin realizar primero pruebas de su rendimiento con un set de testeo, además de que mejoraría en gran medida si se entrenara con un set más equilibrado.

Modelo alternativo

Vistos los resultados obtenidos con el modelo anterior, nos preguntamos si podríamos mejorarlo más. Y es que, como explicamos, el predictor lcavol tiene una baja significancia, esto junto con la alta correlación con el resto de predictores del modelo nos lleva a pensar que podemos mejorar el modelo al eliminar este predictor.

Vamos a proceder de forma similar al modelo original:

modelo_logistico.nolcavol <- glm(svi ~ lcp + lpsa, data = prostate, family = binomial)
# Creamos una lista con los coeficientes preparados para el display
coefs = sapply(modelo_logistico.nolcavol$coefficients, function(x) round(x,2))
summary(modelo_logistico.nolcavol)

## 
## Call:
## glm(formula = svi ~ lcp + lpsa, family = binomial, data = prostate)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -8.5322     2.4718  -3.452 0.000557 ***
## lcp           1.3434     0.3830   3.507 0.000452 ***
## lpsa          2.1982     0.7416   2.964 0.003035 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 101.353  on 96  degrees of freedom
## Residual deviance:  39.653  on 94  degrees of freedom
## AIC: 45.653
## 
## Number of Fisher Scoring iterations: 7

El modelo resulta en la siguiente función: \[ \text{svi} = -8.53 + 1.34 \cdot \text{lcp}+ 2.2 \cdot \text{lpsa} \]

Además, los \(p-\)values nos muestran que todas las variables utilizadas, lcp y lpsa en este caso, son significantes. Vamos a comparar los resultados obtenidos

Vamos a realizar una predicción de nuestro set original usando el modelo:

threshold = 0.5
log.probs = predict(modelo_logistico.nolcavol, type = "response")
log.preds = rep(0, length(log.probs))
log.preds[log.probs > threshold] = 1

Veamos la tabla de valores originales:

table(log.preds)

## log.preds
##  0  1 
## 78 19

Y a continuación la tabla de valores obtenidos:

table(prostate$svi)

## 
##  1  2 
## 76 21

Observamos como hemos conseguido predecir un valor más utilizando este modelo, que se desprende de la variable lcavol, llegando así al 98% de precisión. Sin embargo, como también mencionamos con el modelo anterior, debemos ser críticos con este modelo y prestar especialmente atención al posible overfitting y al desbalance del conjunto de datos.

8.3 Predicción con el modelo

Vamos a utilizar nuestro modelo para predecir un valor de svi para valores concretos de los predictores lcavol, lcp y lpsa. Utilizamos los siguientes valores:

• lcavol = 2.8269,
• lcp = 1.843,
• lpsa = 3.285

Primero creamos el data.frame que contenga estos datos:

pred_data = data.frame(
  lcavol = 2.8269,
  lcp = 1.843,
  lpsa = 3.285
)

Realizamos la predicción usando nuestro modelo:

prob <- predict(modelo_logistico, newdata = pred_data, type = "response")

Obtenemos una probabilidad del 77.1% de que svi\(=1\) con estos valores de los predictores.

Predicción usando el modelo alternativo

Vamos a realizar la misma predicción, pero esta vez utilizando el modelo propuesto en el subapartado anterior, que se desprendía del predictor lcavol por su baja significancia estadística. Dado que los datos ya están preparados, podemos calcular la probabilidad con un solo comando (se ignora el valor de lcavol al realizar la predicción):

prob.nolcavol <- predict(modelo_logistico.nolcavol, newdata = pred_data, type = "response")

Obtenemos una probabilidad del 76.2% de que svi\(=1\) con estos valores de los predictores. Podemos observar que la probabilidad es prácticamente la misma que con el modelo anterior.

9. PCA-PCR

9.1 Biplot y summary del modelo

Primero, creamos nuestro modelo de componentes principales PCA:

# Creamos el subset para el PCA
pca.data <- prostate[,c("lcavol", "lweight", "lbph", "lcp")]
# Añadimos el svi a parte, ya que al ser categorico hay que reconvertirlo a numerico
pca.data["svi"] <- as.numeric(prostate$svi)
# Realizamos el PCA
pca.model <- prcomp(pca.data, scale. = TRUE, center = TRUE)
pca.var <- pca.model$sdev ^2
pca.pve <- pca.var / sum(pca.var)

summary del modelo

summary(pca.model)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     1.5341 1.1862 0.7258 0.66771 0.51651
## Proportion of Variance 0.4707 0.2814 0.1054 0.08917 0.05336
## Cumulative Proportion  0.4707 0.7521 0.8575 0.94664 1.00000

Podemos observar en la primera fila del summary de nuestro modelo, Standard deviation, la cantidad de variabilidad que cada componente principal captura de los datos. Como se puede observar, los valores disminuyen con cada componente, siendo la mayor PC1 y la menor PC5.

Desviación estándar de la primera componente principal

La desviación estándar en un análisis de componentes principales representa la cantidad de variabilidad de los datos que es capturada por esta componente, lo cual es una medida clave en el PCA, ya que su objetivo es el de identificar las direcciones, o componentes principales, que mejor capturen la varianza del conjunto de datos.

El primer componente principal PC1 es, por definición, aquel que caputra la mayor varianza posible, maximizando la varianza en las observaciones proyectadas en esa dirección, con lo cual una gran parte de las diferencias entre observaciones se ven reflejadas a lo largo de esta dirección.

Por tanto, se puede decir que la mayor desviación estándar de esta primera componente es consecuencia del objetivo de un PCA.

Componentes principales

Las siguientes filas del summary ahondan en esta descripción, proporcionándonos valores porcentuales y acumulativos de la descripción de la varianza explicada por cada componente. Representemos estos valores gráficamente:

par(mfrow = c(1, 2))
plot(pca.pve, xlab = "Principal Component",
  ylab = "Proportion of Variance Explained", ylim = c(0, 1),
  type = "b")
plot(cumsum(pca.pve), xlab = "Principal Component",
  ylab = "Cumulative Proportion of Variance Explained",
  ylim = c(0, 1), type = "b")

Como podemos observar, un \(75%\) de la varianza de los datos originales se podrían explicar usando tan solo las dos primeras componentes principales, lo cual sugiere que estas dos componentes podrían ser suficientes para representar los datos, con una pequeña pérdida de información sobre la variación de los datos, llegando al \(85%\) en caso de utilizar tres componentes.

biplot con las dos primeras componentes principales

Realicemos un biplot para ver cómo se expresarían nuestros datos en función de esas dos componentes.

biplot(pca.model, scale = 0)

En este biplot podemos observar los puntos que representan a los pacientes en el espacio definido por las dos primeras componentes principales (PC1 y PC2).

Podemos observar diferentes agrupaciones de puntos, lo que indicaría que estas observaciones tendrían valores similares también en las variables originales.

Además, gracias a las flechas que indican la composición de las componentes principales según las variables originales, podemos observar que, al haber un ángulo pequeño entre ellas, las variables lcavol, lcp y svi están fuertemente correlacionadas, como se corrobora en apartados anteriores. Lo mismo ocurriría con lbph y lweight aunque en menor medida, al haber un ángulo de separación mayor entre ellas.

9.2 PCR

Iniciamos con una carga de las librerias necesarias

library(caret)

Análisis exploratorio de datos y correlación

Iniciamos echando una pequeña vista a la correlación existente entre los datos. Se ha eliminado la diagonal para mejorar la visualización.

df <- prostate
df$svi <- as.numeric(df$svi)

corrplot(cor(df[, c("lcavol", "lweight", "lbph", "lcp", "svi")]), 
         type = "upper", 
         tl.cex = 0.7, 
         diag = FALSE,
         addCoef.col = "black") 

Si hubiéramos usado la función findCorrelation, no se vería ninguna hasta modificar el valor por defecto de cutoff = 0.95.

findCorrelation(cor(df[, c("lcavol", "lweight", "lbph", "lcp", "svi")]), 
                names = TRUE,
                cutoff = 0.44)

## [1] "lcavol"  "lcp"     "lweight"

Ajuste de la Regresión con Componentes Principales (PCR)

suppressPackageStartupMessages({
  library(pls)
})

Para abordar la colinealidad, se ajusta un modelo de PCR y uso validación cruzada para evaluar el error de predicción.

Se incluye scale = TRUE para normalizar los predictores y validation = "CV" para usar validación cruzada en la evaluación:

set.seed(42)
pcr.fit = pcr(lpsa ~ lcavol + lweight + lbph + lcp + svi, 
              data = df, 
              scale = TRUE, 
              validation = "CV")

Visualización del MSE y selección de componentes principales

Para decidir el número óptimo de componentes principales, se traza las puntuaciones de la validación cruzada utilizando validationplot(). Visualmente ya se ve que el menor error de validación cruzada es con M = 5.

Además se ve que el aumento de componentes favorece en los resultados. Haciendo uso completo de los pedidos en el enunciado. Aunque de forma relativa el uso de 1 a 4 componentes podría ser equivalente si quisiéramos usar menos.

validationplot(pcr.fit, val.type = "MSEP")

Fijándonos en el summary() podemos ver mismamente que en M = 5 se captura el 100% de la varianza en los predictores.

summary(pcr.fit)

## Data:    X dimension: 97 5 
##  Y dimension: 97 1
## Fit method: svdpc
## Number of components considered: 5
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps
## CV            1.16   0.7694   0.7558   0.7551   0.7567   0.7281
## adjCV         1.16   0.7679   0.7547   0.7541   0.7548   0.7258
## 
## TRAINING: % variance explained
##       1 comps  2 comps  3 comps  4 comps  5 comps
## X       47.07    75.21    85.75    94.66   100.00
## lpsa    56.48    58.88    59.44    60.56    64.43

O a través de la función which.min() vemos el error más pequeño:

which.min( MSEP( pcr.fit )$val[1,1, ] ) - 1

## 5 comps 
##       5

Comparación con Regresión lineal múltiple

# Predicciones
pcr.pred_full <- predict(pcr.fit, df, ncomp = 5)
multi_pred_full <- predict(modelo_lpsa, newdata = df)



# SSE para PCR
sse_pcr_full <- sum((pcr.pred_full - df$lpsa)^2)

# SSE para Regresión Lineal Múltiple
sse_multi_full <- sum((multi_pred_full - df$lpsa)^2)


n <- nrow(df)
p_pcr <- 5  
p_multi <- length(coef(modelo_lpsa)) 


# Calcular el RSE para PCR
rse_pcr_full <- sqrt(sse_pcr_full / (n - p_pcr))
# Calcular el RSE para Regresión Lineal Múltiple
rse_multi_full <- sqrt(sse_multi_full / (n - p_multi))




# R2 pcr
ss_total_full <- sum((df$lpsa - mean(df$lpsa))^2)
ss_residual_pcr_full <- sum((pcr.pred_full - df$lpsa)^2)
r_squared_pcr_full <- 1 - (ss_residual_pcr_full / ss_total_full)
adj_r_squared_pcr_full <- 1 - ((1 - r_squared_pcr_full) * (nrow(df) - 1) / (nrow(df) - 2 - 1))
# R2 multi
adj_r_squared_multi_full <- summary(modelo_lpsa)$adj.r.squared


cat("PCR: regresión de componentes principales\n",
    "MULTI: Regresión Lineal Múltiple\n\n",
    "############ RSE ############\n\n",
    "RSE (PCR):", rse_pcr_full, "\n",
    "RSE (MULTI):", rse_multi_full, "\n\n",
    "############ R² ############\n\n",
    "R² ajustado (PCR):", adj_r_squared_pcr_full, "\n",
    "R² ajustado (MULTI):", adj_r_squared_multi_full, "\n")

## PCR: regresión de componentes principales
##  MULTI: Regresión Lineal Múltiple
## 
##  ############ RSE ############
## 
##  RSE (PCR): 0.7032094 
##  RSE (MULTI): 0.7070626 
## 
##  ############ R² ############
## 
##  R² ajustado (PCR): 0.6367798 
##  R² ajustado (MULTI): 0.6248055

RSE (Error Estándar de la Regresión):

Ambos modelos, PCR y Regresión Lineal Múltiple (MULTI), tienen un RSE muy similar:

• RSE (PCR): 0.7032
  
• RSE (MULTI): 0.7071

Esto indica que ambos modelos tienen un error estándar de la regresión comparable, lo que sugiere que la dispersión de los errores de las predicciones es casi igual en ambos casos.

R² ajustado:

• El modelo de PCR tiene un R² ajustado de 0.6368, lo que significa que explica el 63.68% de la variabilidad de los datos.
• El modelo de Regresión Lineal Múltiple (MULTI) tiene un R² ajustado de 0.6248, lo que implica que explica el 62.48% de la variabilidad de los datos.

Conclusión:

Aunque ambos modelos presentan RSE similares, lo que refleja una dispersión de errores comparable, el modelo PCR tiene una ligera ventaja en cuanto a R² ajustado, indicando que PCR explica mejor la variabilidad de los datos. Este pequeño beneficio en el ajuste hace que PCR sea el modelo preferido en este caso, ya que demuestra un mayor rendimiento al ajustar los datos.

Sin embargo, la diferencia no es tan significativa, por lo que dependiendo del contexto y los requisitos del análisis, Regresión Lineal Múltiple (MULTI) también puede ser una opción válida si se busca simplicidad o interpretabilidad adicional.